%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the skeleton of C given the outer contour D
%
% inputs:
% kh         Helmholtz parameter.
% C          contour
% matA       discretization matrix corresponding to contour C alone
% acc        accuracy for builsing low rank approximations
% D          outer contour
% h_D        point spacing on D
%
% outputs:
% Is         indices of points in the skeleton
% Cs         subset of C with only skeleton points
% V          projection matrix
% matP       proxy matrix
% D          outer contour (TO REMOVE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Is,Cs,V,matP,D] = get_skeleton(kh,C,matA,acc,D,h_D)

n = size(C,2);
M = size(D,2);

% if kh*h_D/(2*pi) > 1./5
%     error('mesh too coarse')
% end
% fprintf(1,'  h / lambda = %g\n', kh*h_D/(2*pi))

matB = get_B(kh,C,1:n,D); matB = matB.';
matC = get_C(kh,C,1:n,D,h_D);

[T, I] = id_decomp([matB; matC'],acc,'MAT');

ns = size(T,1);
fprintf(1,'Skeleton : %6d -> %6d points. \n',n,ns)

% fprintf(1,'\nInitial # of points      : %d \n',n)
% fprintf(1,'# of points in the skeleton: %d \n',ns)

V = [eye(ns);T'];
V(I,:) = V;

Is = I(1:ns);
Cs = C(:,Is);

matBs = matB(:,Is);
matCs = matC(Is,:);

matP = V' * inv(matA) * V;
% 
% matS  = matB *inv(matA)*matC;
% matSs = matBs*matP     *matCs;
% 
% err = matS - matSs;
% fprintf(1,'|S - Ss|   = %g\n',norm(err,Inf))

return